home *** CD-ROM | disk | FTP | other *** search
/ Mac-Source 1994 July / Mac-Source_July_1994.iso / C and C++ / Science⁄Math / meschach / oldmeschach.shar < prev    next >
Encoding:
Text File  |  1994-06-07  |  21.5 KB  |  870 lines  |  [TEXT/ttxt]
  1. # to unbundle, sh this file (in an empty directory)
  2. echo conjgrad.c 1>&2
  3. sed >conjgrad.c <<'//GO.SYSIN DD conjgrad.c' 's/^-//'
  4. -
  5. -/**************************************************************************
  6. -**
  7. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  8. -**
  9. -**                 Meschach Library
  10. -** 
  11. -** This Meschach Library is provided "as is" without any express 
  12. -** or implied warranty of any kind with respect to this software. 
  13. -** In particular the authors shall not be liable for any direct, 
  14. -** indirect, special, incidental or consequential damages arising 
  15. -** in any way from use of the software.
  16. -** 
  17. -** Everyone is granted permission to copy, modify and redistribute this
  18. -** Meschach Library, provided:
  19. -**  1.  All copies contain this copyright notice.
  20. -**  2.  All modified copies shall carry a notice stating who
  21. -**      made the last modification and the date of such modification.
  22. -**  3.  No charge is made for this software or works derived from it.  
  23. -**      This clause shall not be construed as constraining other software
  24. -**      distributed on the same medium as this software, nor is a
  25. -**      distribution fee considered a charge.
  26. -**
  27. -***************************************************************************/
  28. -
  29. -
  30. -/*
  31. -    Conjugate gradient routines file
  32. -    Uses sparse matrix input & sparse Cholesky factorisation in pccg().
  33. -
  34. -    All the following routines use routines to define a matrix
  35. -        rather than use any explicit representation
  36. -        (with the exeception of the pccg() pre-conditioner)
  37. -    The matrix A is defined by
  38. -
  39. -        VEC *(*A)(void *params, VEC *x, VEC *y)
  40. -
  41. -    where y = A.x on exit, and y is returned. The params argument is
  42. -    intended to make it easier to re-use & modify such routines.
  43. -
  44. -    If we have a sparse matrix data structure
  45. -        SPMAT    *A_mat;
  46. -    then these can be used by passing sp_mv_mlt as the function, and
  47. -    A_mat as the param.
  48. -*/
  49. -
  50. -#include    <stdio.h>
  51. -#include    <math.h>
  52. -#include    "matrix.h"
  53. -#include    "sparse.h"
  54. -static char    rcsid[] = "$Id: conjgrad.c,v 1.4 1994/01/13 05:36:45 des Exp $";
  55. -
  56. -
  57. -/* #define    MAX_ITER    10000 */
  58. -static    int    max_iter = 10000;
  59. -int    cg_num_iters;
  60. -
  61. -/* matrix-as-routine type definition */
  62. -/* #ifdef ANSI_C */
  63. -/* typedef VEC    *(*MTX_FN)(void *params, VEC *x, VEC *out); */
  64. -/* #else */
  65. -typedef VEC    *(*MTX_FN)();
  66. -/* #endif */
  67. -#ifdef ANSI_C
  68. -VEC    *spCHsolve(SPMAT *,VEC *,VEC *);
  69. -#else
  70. -VEC    *spCHsolve();
  71. -#endif
  72. -
  73. -/* cg_set_maxiter -- sets maximum number of iterations if numiter > 1
  74. -    -- just returns current max_iter otherwise
  75. -    -- returns old maximum */
  76. -int    cg_set_maxiter(numiter)
  77. -int    numiter;
  78. -{
  79. -    int    temp;
  80. -
  81. -    if ( numiter < 2 )
  82. -        return max_iter;
  83. -    temp = max_iter;
  84. -    max_iter = numiter;
  85. -    return temp;
  86. -}
  87. -
  88. -
  89. -/* pccg -- solves A.x = b using pre-conditioner M
  90. -            (assumed factored a la spCHfctr())
  91. -    -- results are stored in x (if x != NULL), which is returned */
  92. -VEC    *pccg(A,A_params,M_inv,M_params,b,eps,x)
  93. -MTX_FN    A, M_inv;
  94. -VEC    *b, *x;
  95. -double    eps;
  96. -void    *A_params, *M_params;
  97. -{
  98. -    VEC    *r = VNULL, *p = VNULL, *q = VNULL, *z = VNULL;
  99. -    int    k;
  100. -    Real    alpha, beta, ip, old_ip, norm_b;
  101. -
  102. -    if ( ! A || ! b )
  103. -        error(E_NULL,"pccg");
  104. -    if ( x == b )
  105. -        error(E_INSITU,"pccg");
  106. -    x = v_resize(x,b->dim);
  107. -    if ( eps <= 0.0 )
  108. -        eps = MACHEPS;
  109. -
  110. -    r = v_get(b->dim);
  111. -    p = v_get(b->dim);
  112. -    q = v_get(b->dim);
  113. -    z = v_get(b->dim);
  114. -
  115. -    norm_b = v_norm2(b);
  116. -
  117. -    v_zero(x);
  118. -    r = v_copy(b,r);
  119. -    old_ip = 0.0;
  120. -    for ( k = 0; ; k++ )
  121. -    {
  122. -        if ( v_norm2(r) < eps*norm_b )
  123. -            break;
  124. -        if ( k > max_iter )
  125. -            error(E_ITER,"pccg");
  126. -        if ( M_inv )
  127. -            (*M_inv)(M_params,r,z);
  128. -        else
  129. -            v_copy(r,z);    /* M == identity */
  130. -        ip = in_prod(z,r);
  131. -        if ( k )    /* if ( k > 0 ) ... */
  132. -        {
  133. -            beta = ip/old_ip;
  134. -            p = v_mltadd(z,p,beta,p);
  135. -        }
  136. -        else        /* if ( k == 0 ) ... */
  137. -        {
  138. -            beta = 0.0;
  139. -            p = v_copy(z,p);
  140. -            old_ip = 0.0;
  141. -        }
  142. -        q = (*A)(A_params,p,q);
  143. -        alpha = ip/in_prod(p,q);
  144. -        x = v_mltadd(x,p,alpha,x);
  145. -        r = v_mltadd(r,q,-alpha,r);
  146. -        old_ip = ip;
  147. -    }
  148. -    cg_num_iters = k;
  149. -
  150. -    V_FREE(p);
  151. -    V_FREE(q);
  152. -    V_FREE(r);
  153. -    V_FREE(z);
  154. -
  155. -    return x;
  156. -}
  157. -
  158. -/* sp_pccg -- a simple interface to pccg() which uses sparse matrix
  159. -        data structures
  160. -    -- assumes that LLT contains the Cholesky factorisation of the
  161. -        actual pre-conditioner */
  162. -VEC    *sp_pccg(A,LLT,b,eps,x)
  163. -SPMAT    *A, *LLT;
  164. -VEC    *b, *x;
  165. -double    eps;
  166. -{    return pccg(sp_mv_mlt,A,spCHsolve,LLT,b,eps,x);        }
  167. -
  168. -
  169. -/*
  170. -    Routines for performing the CGS (Conjugate Gradient Squared)
  171. -    algorithm of P. Sonneveld:
  172. -        "CGS, a fast Lanczos-type solver for nonsymmetric linear
  173. -        systems", SIAM J. Sci. & Stat. Comp. v. 10, pp. 36--52
  174. -*/
  175. -
  176. -/* cgs -- uses CGS to compute a solution x to A.x=b
  177. -    -- the matrix A is not passed explicitly, rather a routine
  178. -        A is passed where A(x,Ax,params) computes
  179. -        Ax = A.x
  180. -    -- the computed solution is passed */
  181. -VEC    *cgs(A,A_params,b,r0,tol,x)
  182. -MTX_FN    A;
  183. -VEC    *x, *b;
  184. -VEC    *r0;        /* tilde r0 parameter -- should be random??? */
  185. -double    tol;        /* error tolerance used */
  186. -void    *A_params;
  187. -{
  188. -    VEC    *p, *q, *r, *u, *v, *tmp1, *tmp2;
  189. -    Real    alpha, beta, norm_b, rho, old_rho, sigma;
  190. -    int    iter;
  191. -
  192. -    if ( ! A || ! x || ! b || ! r0 )
  193. -        error(E_NULL,"cgs");
  194. -    if ( x->dim != b->dim || r0->dim != x->dim )
  195. -        error(E_SIZES,"cgs");
  196. -    if ( tol <= 0.0 )
  197. -        tol = MACHEPS;
  198. -
  199. -    p = v_get(x->dim);
  200. -    q = v_get(x->dim);
  201. -    r = v_get(x->dim);
  202. -    u = v_get(x->dim);
  203. -    v = v_get(x->dim);
  204. -    tmp1 = v_get(x->dim);
  205. -    tmp2 = v_get(x->dim);
  206. -
  207. -    norm_b = v_norm2(b);
  208. -    (*A)(A_params,x,tmp1);
  209. -    v_sub(b,tmp1,r);
  210. -    v_zero(p);    v_zero(q);
  211. -    old_rho = 1.0;
  212. -
  213. -    iter = 0;
  214. -    while ( v_norm2(r) > tol*norm_b )
  215. -    {
  216. -        if ( ++iter > max_iter ) break;
  217. -        /*    error(E_ITER,"cgs");  */
  218. -        rho = in_prod(r0,r);
  219. -        if ( old_rho == 0.0 )
  220. -            error(E_SING,"cgs");
  221. -        beta = rho/old_rho;
  222. -        v_mltadd(r,q,beta,u);
  223. -        v_mltadd(q,p,beta,tmp1);
  224. -        v_mltadd(u,tmp1,beta,p);
  225. -
  226. -        (*A)(A_params,p,v);
  227. -
  228. -        sigma = in_prod(r0,v);
  229. -        if ( sigma == 0.0 )
  230. -            error(E_SING,"cgs");
  231. -        alpha = rho/sigma;
  232. -        v_mltadd(u,v,-alpha,q);
  233. -        v_add(u,q,tmp1);
  234. -
  235. -        (*A)(A_params,tmp1,tmp2);
  236. -
  237. -        v_mltadd(r,tmp2,-alpha,r);
  238. -        v_mltadd(x,tmp1,alpha,x);
  239. -
  240. -        old_rho = rho;
  241. -    }
  242. -    cg_num_iters = iter;
  243. -
  244. -    V_FREE(p);    V_FREE(q);    V_FREE(r);
  245. -    V_FREE(u);    V_FREE(v);
  246. -    V_FREE(tmp1);    V_FREE(tmp2);
  247. -
  248. -    return x;
  249. -}
  250. -
  251. -/* sp_cgs -- simple interface for SPMAT data structures */
  252. -VEC    *sp_cgs(A,b,r0,tol,x)
  253. -SPMAT    *A;
  254. -VEC    *b, *r0, *x;
  255. -double    tol;
  256. -{    return cgs(sp_mv_mlt,A,b,r0,tol,x);    }
  257. -
  258. -/*
  259. -    Routine for performing LSQR -- the least squares QR algorithm
  260. -    of Paige and Saunders:
  261. -        "LSQR: an algorithm for sparse linear equations and
  262. -        sparse least squares", ACM Trans. Math. Soft., v. 8
  263. -        pp. 43--71 (1982)
  264. -*/
  265. -/* lsqr -- sparse CG-like least squares routine:
  266. -    -- finds min_x ||A.x-b||_2 using A defined through A & AT
  267. -    -- returns x (if x != NULL) */
  268. -VEC    *lsqr(A,AT,A_params,b,tol,x)
  269. -MTX_FN    A, AT;    /* AT is A transposed */
  270. -VEC    *x, *b;
  271. -double    tol;        /* error tolerance used */
  272. -void    *A_params;
  273. -{
  274. -    VEC    *u, *v, *w, *tmp;
  275. -    Real    alpha, beta, norm_b, phi, phi_bar,
  276. -                rho, rho_bar, rho_max, theta;
  277. -    Real    s, c;    /* for Givens' rotations */
  278. -    int    iter, m, n;
  279. -
  280. -    if ( ! b || ! x )
  281. -        error(E_NULL,"lsqr");
  282. -    if ( tol <= 0.0 )
  283. -        tol = MACHEPS;
  284. -
  285. -    m = b->dim;    n = x->dim;
  286. -    u = v_get((u_int)m);
  287. -    v = v_get((u_int)n);
  288. -    w = v_get((u_int)n);
  289. -    tmp = v_get((u_int)n);
  290. -    norm_b = v_norm2(b);
  291. -
  292. -    v_zero(x);
  293. -    beta = v_norm2(b);
  294. -    if ( beta == 0.0 )
  295. -        return x;
  296. -    sv_mlt(1.0/beta,b,u);
  297. -    tracecatch((*AT)(A_params,u,v),"lsqr");
  298. -    alpha = v_norm2(v);
  299. -    if ( alpha == 0.0 )
  300. -        return x;
  301. -    sv_mlt(1.0/alpha,v,v);
  302. -    v_copy(v,w);
  303. -    phi_bar = beta;        rho_bar = alpha;
  304. -
  305. -    rho_max = 1.0;
  306. -    iter = 0;
  307. -    do {
  308. -        if ( ++iter > max_iter )
  309. -            error(E_ITER,"lsqr");
  310. -
  311. -        tmp = v_resize(tmp,m);
  312. -        tracecatch((*A) (A_params,v,tmp),"lsqr");
  313. -
  314. -        v_mltadd(tmp,u,-alpha,u);
  315. -        beta = v_norm2(u);    sv_mlt(1.0/beta,u,u);
  316. -
  317. -        tmp = v_resize(tmp,n);
  318. -        tracecatch((*AT)(A_params,u,tmp),"lsqr");
  319. -        v_mltadd(tmp,v,-beta,v);
  320. -        alpha = v_norm2(v);    sv_mlt(1.0/alpha,v,v);
  321. -
  322. -        rho = sqrt(rho_bar*rho_bar+beta*beta);
  323. -        if ( rho > rho_max )
  324. -            rho_max = rho;
  325. -        c   = rho_bar/rho;
  326. -        s   = beta/rho;
  327. -        theta   =  s*alpha;
  328. -        rho_bar = -c*alpha;
  329. -        phi     =  c*phi_bar;
  330. -        phi_bar =  s*phi_bar;
  331. -
  332. -        /* update x & w */
  333. -        if ( rho == 0.0 )
  334. -            error(E_SING,"lsqr");
  335. -        v_mltadd(x,w,phi/rho,x);
  336. -        v_mltadd(v,w,-theta/rho,w);
  337. -    } while ( fabs(phi_bar*alpha*c) > tol*norm_b/rho_max );
  338. -
  339. -    cg_num_iters = iter;
  340. -
  341. -    V_FREE(tmp);    V_FREE(u);    V_FREE(v);    V_FREE(w);
  342. -
  343. -    return x;
  344. -}
  345. -
  346. -/* sp_lsqr -- simple interface for SPMAT data structures */
  347. -VEC    *sp_lsqr(A,b,tol,x)
  348. -SPMAT    *A;
  349. -VEC    *b, *x;
  350. -double    tol;
  351. -{    return lsqr(sp_mv_mlt,sp_vm_mlt,A,b,tol,x);    }
  352. -
  353. //GO.SYSIN DD conjgrad.c
  354. echo lanczos.c 1>&2
  355. sed >lanczos.c <<'//GO.SYSIN DD lanczos.c' 's/^-//'
  356. -
  357. -/**************************************************************************
  358. -**
  359. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  360. -**
  361. -**                 Meschach Library
  362. -** 
  363. -** This Meschach Library is provided "as is" without any express 
  364. -** or implied warranty of any kind with respect to this software. 
  365. -** In particular the authors shall not be liable for any direct, 
  366. -** indirect, special, incidental or consequential damages arising 
  367. -** in any way from use of the software.
  368. -** 
  369. -** Everyone is granted permission to copy, modify and redistribute this
  370. -** Meschach Library, provided:
  371. -**  1.  All copies contain this copyright notice.
  372. -**  2.  All modified copies shall carry a notice stating who
  373. -**      made the last modification and the date of such modification.
  374. -**  3.  No charge is made for this software or works derived from it.  
  375. -**      This clause shall not be construed as constraining other software
  376. -**      distributed on the same medium as this software, nor is a
  377. -**      distribution fee considered a charge.
  378. -**
  379. -***************************************************************************/
  380. -
  381. -
  382. -/*
  383. -    File containing Lanczos type routines for finding eigenvalues
  384. -    of large, sparse, symmetic matrices
  385. -*/
  386. -
  387. -#include    <stdio.h>
  388. -#include    <math.h>
  389. -#include    "matrix.h"
  390. -#include    "sparse.h"
  391. -
  392. -static char rcsid[] = "$Id: lanczos.c,v 1.4 1994/01/13 05:28:24 des Exp $";
  393. -
  394. -#ifdef ANSI_C
  395. -extern    VEC    *trieig(VEC *,VEC *,MAT *);
  396. -#else
  397. -extern    VEC    *trieig();
  398. -#endif
  399. -
  400. -/* lanczos -- raw lanczos algorithm -- no re-orthogonalisation
  401. -    -- creates T matrix of size == m,
  402. -        but no larger than before beta_k == 0
  403. -    -- uses passed routine to do matrix-vector multiplies */
  404. -void    lanczos(A_fn,A_params,m,x0,a,b,beta2,Q)
  405. -VEC    *(*A_fn)();    /* VEC *(*A_fn)(void *A_params,VEC *in, VEC *out) */
  406. -void    *A_params;
  407. -int    m;
  408. -VEC    *x0, *a, *b;
  409. -Real    *beta2;
  410. -MAT    *Q;
  411. -{
  412. -    int    j;
  413. -    VEC    *v, *w, *tmp;
  414. -    Real    alpha, beta;
  415. -
  416. -    if ( ! A_fn || ! x0 || ! a || ! b )
  417. -        error(E_NULL,"lanczos");
  418. -    if ( m <= 0 )
  419. -        error(E_BOUNDS,"lanczos");
  420. -    if ( Q && ( Q->m < x0->dim || Q->n < m ) )
  421. -        error(E_SIZES,"lanczos");
  422. -
  423. -    a = v_resize(a,(u_int)m);    b = v_resize(b,(u_int)(m-1));
  424. -    v = v_get(x0->dim);
  425. -    w = v_get(x0->dim);
  426. -    tmp = v_get(x0->dim);
  427. -
  428. -    beta = 1.0;
  429. -    /* normalise x0 as w */
  430. -    sv_mlt(1.0/v_norm2(x0),x0,w);
  431. -
  432. -    (*A_fn)(A_params,w,v);
  433. -
  434. -    for ( j = 0; j < m; j++ )
  435. -    {
  436. -        /* store w in Q if Q not NULL */
  437. -        if ( Q )
  438. -            set_col(Q,j,w);
  439. -
  440. -        alpha = in_prod(w,v);
  441. -        a->ve[j] = alpha;
  442. -        v_mltadd(v,w,-alpha,v);
  443. -        beta = v_norm2(v);
  444. -        if ( beta == 0.0 )
  445. -        {
  446. -            v_resize(a,(u_int)j+1);
  447. -            v_resize(b,(u_int)j);
  448. -            *beta2 = 0.0;
  449. -            if ( Q )
  450. -            Q = m_resize(Q,Q->m,j+1);
  451. -            return;
  452. -        }
  453. -        if ( j < m-1 )
  454. -            b->ve[j] = beta;
  455. -        v_copy(w,tmp);
  456. -        sv_mlt(1/beta,v,w);
  457. -        sv_mlt(-beta,tmp,v);
  458. -        (*A_fn)(A_params,w,tmp);
  459. -        v_add(v,tmp,v);
  460. -    }
  461. -    *beta2 = beta;
  462. -
  463. -
  464. -    V_FREE(v);    V_FREE(w);    V_FREE(tmp);
  465. -}
  466. -
  467. -extern    double    frexp(), ldexp();
  468. -
  469. -/* product -- returns the product of a long list of numbers
  470. -    -- answer stored in mant (mantissa) and expt (exponent) */
  471. -static    double    product(a,offset,expt)
  472. -VEC    *a;
  473. -double    offset;
  474. -int    *expt;
  475. -{
  476. -    Real    mant, tmp_fctr;
  477. -    int    i, tmp_expt;
  478. -
  479. -    if ( ! a )
  480. -        error(E_NULL,"product");
  481. -
  482. -    mant = 1.0;
  483. -    *expt = 0;
  484. -    if ( offset == 0.0 )
  485. -        for ( i = 0; i < a->dim; i++ )
  486. -        {
  487. -            mant *= frexp(a->ve[i],&tmp_expt);
  488. -            *expt += tmp_expt;
  489. -            if ( ! (i % 10) )
  490. -            {
  491. -                mant = frexp(mant,&tmp_expt);
  492. -                *expt += tmp_expt;
  493. -            }
  494. -        }
  495. -    else
  496. -        for ( i = 0; i < a->dim; i++ )
  497. -        {
  498. -            tmp_fctr = a->ve[i] - offset;
  499. -            tmp_fctr += (tmp_fctr > 0.0 ) ? -MACHEPS*offset :
  500. -                             MACHEPS*offset;
  501. -            mant *= frexp(tmp_fctr,&tmp_expt);
  502. -            *expt += tmp_expt;
  503. -            if ( ! (i % 10) )
  504. -            {
  505. -                mant = frexp(mant,&tmp_expt);
  506. -                *expt += tmp_expt;
  507. -            }
  508. -        }
  509. -
  510. -    mant = frexp(mant,&tmp_expt);
  511. -    *expt += tmp_expt;
  512. -
  513. -    return mant;
  514. -}
  515. -
  516. -/* product2 -- returns the product of a long list of numbers
  517. -    -- answer stored in mant (mantissa) and expt (exponent) */
  518. -static    double    product2(a,k,expt)
  519. -VEC    *a;
  520. -int    k;    /* entry of a to leave out */
  521. -int    *expt;
  522. -{
  523. -    Real    mant, mu, tmp_fctr;
  524. -    int    i, tmp_expt;
  525. -
  526. -    if ( ! a )
  527. -        error(E_NULL,"product2");
  528. -    if ( k < 0 || k >= a->dim )
  529. -        error(E_BOUNDS,"product2");
  530. -
  531. -    mant = 1.0;
  532. -    *expt = 0;
  533. -    mu = a->ve[k];
  534. -    for ( i = 0; i < a->dim; i++ )
  535. -    {
  536. -        if ( i == k )
  537. -            continue;
  538. -        tmp_fctr = a->ve[i] - mu;
  539. -        tmp_fctr += ( tmp_fctr > 0.0 ) ? -MACHEPS*mu : MACHEPS*mu;
  540. -        mant *= frexp(tmp_fctr,&tmp_expt);
  541. -        *expt += tmp_expt;
  542. -        if ( ! (i % 10) )
  543. -        {
  544. -            mant = frexp(mant,&tmp_expt);
  545. -            *expt += tmp_expt;
  546. -        }
  547. -    }
  548. -    mant = frexp(mant,&tmp_expt);
  549. -    *expt += tmp_expt;
  550. -
  551. -    return mant;
  552. -}
  553. -
  554. -/* dbl_cmp -- comparison function to pass to qsort() */
  555. -static    int    dbl_cmp(x,y)
  556. -Real    *x, *y;
  557. -{
  558. -    Real    tmp;
  559. -
  560. -    tmp = *x - *y;
  561. -    return (tmp > 0 ? 1 : tmp < 0 ? -1: 0);
  562. -}
  563. -
  564. -/* lanczos2 -- lanczos + error estimate for every e-val
  565. -    -- uses Cullum & Willoughby approach, Sparse Matrix Proc. 1978
  566. -    -- returns multiple e-vals where multiple e-vals may not exist
  567. -    -- returns evals vector */
  568. -VEC    *lanczos2(A_fn,A_params,m,x0,evals,err_est)
  569. -VEC    *(*A_fn)();
  570. -void    *A_params;
  571. -int    m;
  572. -VEC    *x0;        /* initial vector */
  573. -VEC    *evals;        /* eigenvalue vector */
  574. -VEC    *err_est;    /* error estimates of eigenvalues */
  575. -{
  576. -    VEC        *a;
  577. -    static    VEC    *b=VNULL, *a2=VNULL, *b2=VNULL;
  578. -    Real    beta, pb_mant, det_mant, det_mant1, det_mant2;
  579. -    int    i, pb_expt, det_expt, det_expt1, det_expt2;
  580. -
  581. -    if ( ! A_fn || ! x0 )
  582. -        error(E_NULL,"lanczos2");
  583. -    if ( m <= 0 )
  584. -        error(E_RANGE,"lanczos2");
  585. -
  586. -    a = evals;
  587. -    a = v_resize(a,(u_int)m);
  588. -    b = v_resize(b,(u_int)(m-1));
  589. -    MEM_STAT_REG(b,TYPE_VEC);
  590. -
  591. -    lanczos(A_fn,A_params,m,x0,a,b,&beta,MNULL);
  592. -
  593. -    /* printf("# beta =%g\n",beta); */
  594. -    pb_mant = 0.0;
  595. -    if ( err_est )
  596. -    {
  597. -        pb_mant = product(b,(double)0.0,&pb_expt);
  598. -        /* printf("# pb_mant = %g, pb_expt = %d\n",pb_mant, pb_expt); */
  599. -    }
  600. -
  601. -    /* printf("# diags =\n");    out_vec(a); */
  602. -    /* printf("# off diags =\n");    out_vec(b); */
  603. -    a2 = v_resize(a2,a->dim - 1);
  604. -    b2 = v_resize(b2,b->dim - 1);
  605. -    MEM_STAT_REG(a2,TYPE_VEC);
  606. -    MEM_STAT_REG(b2,TYPE_VEC);
  607. -    for ( i = 0; i < a2->dim - 1; i++ )
  608. -    {
  609. -        a2->ve[i] = a->ve[i+1];
  610. -        b2->ve[i] = b->ve[i+1];
  611. -    }
  612. -    a2->ve[a2->dim-1] = a->ve[a2->dim];
  613. -
  614. -    trieig(a,b,MNULL);
  615. -
  616. -    /* sort evals as a courtesy */
  617. -    qsort((void *)(a->ve),(int)(a->dim),sizeof(Real),(int (*)())dbl_cmp);
  618. -
  619. -    /* error estimates */
  620. -    if ( err_est )
  621. -    {
  622. -        err_est = v_resize(err_est,(u_int)m);
  623. -
  624. -        trieig(a2,b2,MNULL);
  625. -        /* printf("# a =\n");    out_vec(a); */
  626. -        /* printf("# a2 =\n");    out_vec(a2); */
  627. -
  628. -        for ( i = 0; i < a->dim; i++ )
  629. -        {
  630. -            det_mant1 = product2(a,i,&det_expt1);
  631. -            det_mant2 = product(a2,(double)a->ve[i],&det_expt2);
  632. -            /* printf("# det_mant1=%g, det_expt1=%d\n",
  633. -                    det_mant1,det_expt1); */
  634. -            /* printf("# det_mant2=%g, det_expt2=%d\n",
  635. -                    det_mant2,det_expt2); */
  636. -            if ( det_mant1 == 0.0 )
  637. -            {   /* multiple e-val of T */
  638. -                err_est->ve[i] = 0.0;
  639. -                continue;
  640. -            }
  641. -            else if ( det_mant2 == 0.0 )
  642. -            {
  643. -                err_est->ve[i] = HUGE;
  644. -                continue;
  645. -            }
  646. -            if ( (det_expt1 + det_expt2) % 2 )
  647. -                /* if odd... */
  648. -                det_mant = sqrt(2.0*fabs(det_mant1*det_mant2));
  649. -            else /* if even... */
  650. -                det_mant = sqrt(fabs(det_mant1*det_mant2));
  651. -            det_expt = (det_expt1+det_expt2)/2;
  652. -            err_est->ve[i] = fabs(beta*
  653. -                ldexp(pb_mant/det_mant,pb_expt-det_expt));
  654. -        }
  655. -    }
  656. -
  657. -    return a;
  658. -}
  659. -
  660. -/* sp_lanczos -- version that uses sparse matrix data structure */
  661. -void    sp_lanczos(A,m,x0,a,b,beta2,Q)
  662. -SPMAT    *A;
  663. -int     m;
  664. -VEC     *x0, *a, *b;
  665. -Real  *beta2;
  666. -MAT     *Q;
  667. -{    lanczos(sp_mv_mlt,A,m,x0,a,b,beta2,Q);    }
  668. -
  669. -/* sp_lanczos2 -- version of lanczos2() that uses sparse matrix data
  670. -                    structure */
  671. -VEC    *sp_lanczos2(A,m,x0,evals,err_est)
  672. -SPMAT    *A;
  673. -int    m;
  674. -VEC    *x0;        /* initial vector */
  675. -VEC    *evals;        /* eigenvalue vector */
  676. -VEC    *err_est;    /* error estimates of eigenvalues */
  677. -{    return lanczos2(sp_mv_mlt,A,m,x0,evals,err_est);    }
  678. -
  679. //GO.SYSIN DD lanczos.c
  680. echo arnoldi.c 1>&2
  681. sed >arnoldi.c <<'//GO.SYSIN DD arnoldi.c' 's/^-//'
  682. -
  683. -/**************************************************************************
  684. -**
  685. -** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  686. -**
  687. -**                 Meschach Library
  688. -** 
  689. -** This Meschach Library is provided "as is" without any express 
  690. -** or implied warranty of any kind with respect to this software. 
  691. -** In particular the authors shall not be liable for any direct, 
  692. -** indirect, special, incidental or consequential damages arising 
  693. -** in any way from use of the software.
  694. -** 
  695. -** Everyone is granted permission to copy, modify and redistribute this
  696. -** Meschach Library, provided:
  697. -**  1.  All copies contain this copyright notice.
  698. -**  2.  All modified copies shall carry a notice stating who
  699. -**      made the last modification and the date of such modification.
  700. -**  3.  No charge is made for this software or works derived from it.  
  701. -**      This clause shall not be construed as constraining other software
  702. -**      distributed on the same medium as this software, nor is a
  703. -**      distribution fee considered a charge.
  704. -**
  705. -***************************************************************************/
  706. -
  707. -/*
  708. -    Arnoldi method for finding eigenvalues of large non-symmetric
  709. -        matrices
  710. -*/
  711. -#include    <stdio.h>
  712. -#include    <math.h>
  713. -#include    "matrix.h"
  714. -#include    "matrix2.h"
  715. -#include    "sparse.h"
  716. -
  717. -static char rcsid[] = "$Id: arnoldi.c,v 1.3 1994/01/13 05:45:40 des Exp $";
  718. -
  719. -/* arnoldi -- an implementation of the Arnoldi method */
  720. -MAT    *arnoldi(A,A_param,x0,m,h_rem,Q,H)
  721. -VEC    *(*A)();
  722. -void    *A_param;
  723. -VEC    *x0;
  724. -int    m;
  725. -Real    *h_rem;
  726. -MAT    *Q, *H;
  727. -{
  728. -    static VEC    *v=VNULL, *u=VNULL, *r=VNULL, *s=VNULL, *tmp=VNULL;
  729. -    int    i;
  730. -    Real    h_val;
  731. -
  732. -    if ( ! A || ! Q || ! x0 )
  733. -        error(E_NULL,"arnoldi");
  734. -    if ( m <= 0 )
  735. -        error(E_BOUNDS,"arnoldi");
  736. -    if ( Q->n != x0->dim ||    Q->m != m )
  737. -        error(E_SIZES,"arnoldi");
  738. -
  739. -    m_zero(Q);
  740. -    H = m_resize(H,m,m);
  741. -    m_zero(H);
  742. -    u = v_resize(u,x0->dim);
  743. -    v = v_resize(v,x0->dim);
  744. -    r = v_resize(r,m);
  745. -    s = v_resize(s,m);
  746. -    tmp = v_resize(tmp,x0->dim);
  747. -    MEM_STAT_REG(u,TYPE_VEC);
  748. -    MEM_STAT_REG(v,TYPE_VEC);
  749. -    MEM_STAT_REG(r,TYPE_VEC);
  750. -    MEM_STAT_REG(s,TYPE_VEC);
  751. -    MEM_STAT_REG(tmp,TYPE_VEC);
  752. -    sv_mlt(1.0/v_norm2(x0),x0,v);
  753. -    for ( i = 0; i < m; i++ )
  754. -    {
  755. -        set_row(Q,i,v);
  756. -        u = (*A)(A_param,v,u);
  757. -        r = mv_mlt(Q,u,r);
  758. -        tmp = vm_mlt(Q,r,tmp);
  759. -        v_sub(u,tmp,u);
  760. -        h_val = v_norm2(u);
  761. -        /* if u == 0 then we have an exact subspace */
  762. -        if ( h_val == 0.0 )
  763. -        {
  764. -        *h_rem = h_val;
  765. -        return H;
  766. -        }
  767. -        /* iterative refinement -- ensures near orthogonality */
  768. -        do {
  769. -        s = mv_mlt(Q,u,s);
  770. -        tmp = vm_mlt(Q,s,tmp);
  771. -        v_sub(u,tmp,u);
  772. -        v_add(r,s,r);
  773. -        } while ( v_norm2(s) > 0.1*(h_val = v_norm2(u)) );
  774. -        /* now that u is nearly orthogonal to Q, update H */
  775. -        set_col(H,i,r);
  776. -        if ( i == m-1 )
  777. -        {
  778. -        *h_rem = h_val;
  779. -        continue;
  780. -        }
  781. -        /* H->me[i+1][i] = h_val; */
  782. -        m_set_val(H,i+1,i,h_val);
  783. -        sv_mlt(1.0/h_val,u,v);
  784. -    }
  785. -
  786. -    return H;
  787. -}
  788. -
  789. -/* sp_arnoldi -- uses arnoldi() with an explicit representation of A */
  790. -MAT    *sp_arnoldi(A,x0,m,h_rem,Q,H)
  791. -SPMAT    *A;
  792. -VEC    *x0;
  793. -int    m;
  794. -Real    *h_rem;
  795. -MAT    *Q, *H;
  796. -{    return arnoldi(sp_mv_mlt,A,x0,m,h_rem,Q,H);    }
  797. -
  798. -/* gmres -- generalised minimum residual algorithm of Saad & Schultz
  799. -        SIAM J. Sci. Stat. Comp. v.7, pp.856--869 (1986)
  800. -    -- y is overwritten with the solution */
  801. -VEC    *gmres(A,A_param,m,Q,R,b,tol,x)
  802. -VEC    *(*A)();
  803. -void    *A_param;
  804. -VEC    *b, *x;
  805. -int    m;
  806. -MAT    *Q, *R;
  807. -double    tol;
  808. -{
  809. -    static VEC    *v=VNULL, *u=VNULL, *r=VNULL, *tmp=VNULL, *rhs=VNULL;
  810. -    static VEC    *diag=VNULL, *beta=VNULL;
  811. -    int    i;
  812. -    Real    h_val, norm_b;
  813. -    
  814. -    if ( ! A || ! Q || ! b || ! R )
  815. -    error(E_NULL,"gmres");
  816. -    if ( m <= 0 )
  817. -    error(E_BOUNDS,"gmres");
  818. -    if ( Q->n != b->dim || Q->m != m )
  819. -    error(E_SIZES,"gmres");
  820. -    
  821. -    x = v_copy(b,x);
  822. -    m_zero(Q);
  823. -    R = m_resize(R,m+1,m);
  824. -    m_zero(R);
  825. -    u = v_resize(u,x->dim);
  826. -    v = v_resize(v,x->dim);
  827. -    tmp = v_resize(tmp,x->dim);
  828. -    rhs = v_resize(rhs,m+1);
  829. -    MEM_STAT_REG(u,TYPE_VEC);
  830. -    MEM_STAT_REG(v,TYPE_VEC);
  831. -    MEM_STAT_REG(r,TYPE_VEC);
  832. -    MEM_STAT_REG(tmp,TYPE_VEC);
  833. -    MEM_STAT_REG(rhs,TYPE_VEC);
  834. -    norm_b = v_norm2(x);
  835. -    if ( norm_b == 0.0 )
  836. -    error(E_RANGE,"gmres");
  837. -    sv_mlt(1.0/norm_b,x,v);
  838. -    
  839. -    for ( i = 0; i < m; i++ )
  840. -    {
  841. -    set_row(Q,i,v);
  842. -    tracecatch(u = (*A)(A_param,v,u),"gmres");
  843. -    r = mv_mlt(Q,u,r);
  844. -    tmp = vm_mlt(Q,r,tmp);
  845. -    v_sub(u,tmp,u);
  846. -    h_val = v_norm2(u);
  847. -    set_col(R,i,r);
  848. -    R->me[i+1][i] = h_val;
  849. -    sv_mlt(1.0/h_val,u,v);
  850. -    }
  851. -    
  852. -    /* use i x i submatrix of R */
  853. -    R = m_resize(R,i+1,i);
  854. -    rhs = v_resize(rhs,i+1);
  855. -    v_zero(rhs);
  856. -    rhs->ve[0] = norm_b;
  857. -    tmp = v_resize(tmp,i);
  858. -    diag = v_resize(diag,i+1);
  859. -    beta = v_resize(beta,i+1);
  860. -    MEM_STAT_REG(beta,TYPE_VEC);
  861. -    MEM_STAT_REG(diag,TYPE_VEC);
  862. -    QRfactor(R,diag /* ,beta */);
  863. -    tmp = QRsolve(R,diag, /* beta, */ rhs,tmp);
  864. -    v_resize(tmp,m);
  865. -    vm_mlt(Q,tmp,x);
  866. -
  867. -    return x;
  868. -}
  869. //GO.SYSIN DD arnoldi.c
  870.